In [1]:
%pylab
import numpy as np
import pylab as plb
import tifffile
Using matplotlib backend: WXAgg
Populating the interactive namespace from numpy and matplotlib

In [2]:
#these functions are used for the image overlays.
def overlay_hue_saturation(overlay,bkgnd,
                       gama = 1.0,
                       color_map = plb.cm.jet,
                       alpha = 0.5):
    from skimage import data, color, io, img_as_float
    gama_correct = lambda x: x**gama
    if len(np.shape(bkgnd)) == 2:
        img_color = np.dstack((bkgnd, bkgnd, bkgnd))
    else:
        img_color = bkgnd
        
    if len(np.shape(overlay)) == 2:
        overlay_color = color_map(gama_correct(overlay))[:,:,:3]
        mask = overlay>0
    else:
        overlay_color = overlay
        mask = np.sum(overlay.astype(uint64),axis = 2)>0
    img_hsv = color.rgb2hsv(img_color)
    overlay_hsv = color.rgb2hsv(overlay_color)
    img_hsv[..., 0] += mask*overlay_hsv[..., 0] * alpha
    img_hsv[..., 1] += mask*overlay_hsv[..., 1] * alpha
    img_overlay = color.hsv2rgb(img_hsv)
    return(img_overlay)

def colorized_image(input_img,
                    map_point = 0.5,
                    color_map = plb.cm.jet,
                    alpha = 0.2,
                    gama = 1.0,
                    gain = 1.0):
    from skimage import data, color, io, img_as_float
    gama_correct = lambda x: x**gama
    img = gain*input_img
    img = gama_correct(img)
    img[img>255] = 255
    img = img.astype(uint8)
    bg_color = np.dstack((img, img, img))
    
    mask = img>0
    map_img = mask*map_point
    
    overlay_color = color_map(map_img)[:,:,:3]
    img_hsv = color.rgb2hsv(bg_color)
    
    overlay_hsv = color.rgb2hsv(overlay_color)
    #make hue and saturation of img map to color
    img_hsv[..., 0] += mask*overlay_hsv[..., 0] * alpha
    img_hsv[..., 1] += mask*overlay_hsv[..., 1] * alpha
    img_overlay = color.hsv2rgb(img_hsv)
    return(img_overlay)
In [3]:
#load images
bkgnd_img = tifffile.TiffFile('./registration/65G06_brightfield_cuticle.tif').asarray()
b1_stack = tifffile.TiffFile('./registration/combined_stack_b1.tiff').asarray()
b2_stack = tifffile.TiffFile('./registration/combined_stack_b2.tiff').asarray()
b3_stack = tifffile.TiffFile('./registration/combined_stack_b3.tiff').asarray()

i1_stack = tifffile.TiffFile('./registration/combined_stack_i1.tiff').asarray()
i2_stack = tifffile.TiffFile('./registration/combined_stack_i2.tiff').asarray()

iii1_stack = tifffile.TiffFile('./registration/combined_stack_iii1.tiff').asarray()
iii3_stack = tifffile.TiffFile('./registration/combined_stack_iii3.tiff').asarray()
iii24_stack = tifffile.TiffFile('./registration/combined_stack_iii24.tiff').asarray()

hg1_stack = tifffile.TiffFile('./registration/combined_stack_hg1.tiff').asarray()
hg2_stack = tifffile.TiffFile('./registration/combined_stack_hg2.tiff').asarray()
hg3_stack = tifffile.TiffFile('./registration/combined_stack_hg3.tiff').asarray()
hg4_stack = tifffile.TiffFile('./registration/combined_stack_hg4.tiff').asarray()

tpd_stack = tifffile.TiffFile('./registration/combined_stack_tpd.tiff').asarray()
tpv_stack = tifffile.TiffFile('./registration/combined_stack_tpv.tiff').asarray()

pr_stack = tifffile.TiffFile('./registration/combined_stack_pr.tiff').asarray()

#make max projections
b1_img = np.max(b1_stack,axis = 0)
b2_img = np.max(b2_stack,axis = 0)
b3_img = np.max(b3_stack,axis = 0)

i1_img = np.max(i1_stack,axis = 0)
i2_img = np.max(i2_stack,axis = 0)

iii1_img = np.max(iii1_stack,axis = 0)
iii3_img = np.max(iii3_stack,axis = 0)
iii24_img = np.max(iii24_stack,axis = 0)

hg1_img = np.max(hg1_stack,axis = 0)
hg2_img = np.max(hg2_stack,axis = 0)
hg3_img = np.max(hg3_stack,axis = 0)
hg4_img = np.max(hg4_stack,axis = 0)

tpd_img = np.max(tpd_stack,axis = 0)
tpv_img = np.max(tpv_stack,axis = 0)

pr_img = np.max(pr_stack,axis = 0)
In [4]:
#colorize each muscle

b_map = plb.cm.Blues_r
i_map = plb.cm.Purples_r
iii_map = plb.cm.Reds_r
hg_map = plb.cm.Oranges_r
tp_map = plb.cm.Greens_r
pr_map = plb.cm.Greens_r

cvls = {'b1':0.01,'b2':0.2,'b3':0.4,
        'i1':0.01,'i2':0.2,
        'iii1':0.01,'iii3':0.2,'iii24':0.4,
        'hg1':0.01,'hg2':0.2,'hg3':0.4,'hg4':0.6,
        'tpd':0.01,'tpv':0.1,
        'pr':0.01}

b1_colorized = colorized_image(b1_img,color_map = b_map,map_point = cvls['b1'],alpha = 1.0)
b2_colorized = colorized_image(b2_img,color_map = b_map,map_point = cvls['b2'],alpha = 1.0)
b3_colorized = colorized_image(b3_img,color_map = b_map,map_point = cvls['b3'],alpha = 1.0,gama = 0.65,gain = 40.0)

i1_colorized = colorized_image(i1_img,color_map = i_map,map_point = cvls['i1'],alpha = 1.0)
i2_colorized = colorized_image(i2_img,color_map = i_map,map_point = cvls['i2'],alpha = 1.0,gama = 0.7,gain = 10.0)

iii1_colorized = colorized_image(iii1_img,color_map = iii_map,map_point = cvls['iii1'],alpha = 1.0)
iii3_colorized = colorized_image(iii3_img,color_map = iii_map,map_point = cvls['iii3'],alpha = 1.0)
iii24_colorized = colorized_image(iii24_img,color_map = iii_map,map_point = cvls['iii24'],alpha = 1.0)

hg1_colorized = colorized_image(hg1_img,color_map = hg_map,map_point = cvls['hg1'],alpha = 1.0,gain = 2.0)
hg2_colorized = colorized_image(hg2_img,color_map = hg_map,map_point = cvls['hg2'],alpha = 1.0)
hg3_colorized = colorized_image(hg3_img,color_map = hg_map,map_point = cvls['hg3'],alpha = 1.0)
hg4_colorized = colorized_image(hg4_img,color_map = hg_map,map_point = cvls['hg4'],alpha = 1.0)

tpd_colorized = colorized_image(tpd_img,color_map = tp_map,map_point = cvls['tpd'],alpha = 1.0)
tpv_colorized = colorized_image(tpv_img,color_map = tp_map,map_point = cvls['tpv'],alpha = 1.0)

pr_colorized = colorized_image(pr_img,color_map = pr_map,map_point = cvls['pr'],alpha = 1.0)

#combine the images using addition
img_sum = (tpd_colorized+
           tpv_colorized+
           b1_colorized+
           b2_colorized+
           b3_colorized+
           i1_colorized+
           i2_colorized+
           iii1_colorized+
           iii3_colorized+
           iii24_colorized+
           hg1_colorized+
           hg2_colorized+
           hg3_colorized+
           hg4_colorized+
           pr_colorized)*255

img_sum[img_sum>255] = 255
img_sum = img_sum.astype(uint8)
imshow(overlay_hue_saturation(img_sum[:,:,:3],bkgnd_img,alpha = 1.0))
display(gcf());close()

I recursively apply the function 'mask_stack' to a list of colorized muscle stacks to create a visualization of the muscle population. The order of muscles in the list is important since this determines the way the muscles overlap in the final image. The colorized image of the muscles are then mixed into the greyscale image of the cuticle.

In [5]:
def mask_stack(top_img,bottom_img,alpha = 0.5,gain = 15):
    from skimage import data, color, io, img_as_float
    top_hsv = color.rgb2hsv(top_img)
    mask = top_hsv[...,2]/np.max(top_hsv[...,2])
    mask *= gain
    mask[mask>1] = 1
    mask[mask<1] = 0
    mask = np.dstack((mask,mask,mask)) 
    inverse = 1-mask*alpha
    return top_img*mask + bottom_img*inverse

stacked_muscles = reduce(mask_stack, 
                          [tpd_colorized,
                           tpv_colorized,
                           pr_colorized,
                           b2_colorized,
                           i1_colorized,
                           hg4_colorized,
                           iii24_colorized,
                           i2_colorized,
                           hg2_colorized,
                           hg3_colorized,
                           iii3_colorized,
                           iii1_colorized,
                           b1_colorized,
                           b3_colorized,
                           hg1_colorized])
                            
gama = 1.5
gama_correct = lambda x: x**gama
cut_gamma = gama_correct(bkgnd_img)+600
cuticle = (np.dstack((cut_gamma,cut_gamma,cut_gamma))*0.03).astype(uint8)
muscles = (stacked_muscles*0.9*255)
sumimg = cuticle+muscles
sumimg[sumimg>255] = 255
sumimg = sumimg.astype(uint8)
imshow(sumimg)
display(gcf());close()
In [6]:
tifffile.imsave('stacked_muscles.tiff',uint8(stacked_muscles*255))
tifffile.imsave('muscle_viz.tiff',sumimg)

To define the regions of interest I find the 2d contours around each muscle. For simplicity I then resample the contours to 20 points: hopefully this should allow us to average in the map extracted from other confocal stacks in the future. We then save the lists of contours to use in the construction of the model.

In [7]:
mimgs = {'b1':b1_img,
         'b2':b2_img,
         'b3':b3_img,
         'i1':i1_img,
         'i2':i2_img,
         'iii1':iii1_img,
         'iii3':iii3_img,
         'iii4':iii24_img,
         'hg1':hg1_img,
         'hg2':hg2_img,
         'hg3':hg3_img,
         'hg4':hg4_img,
         'tpv':tpv_img,
         'tpd':tpd_img,
         'pr':pr_img}
model_data = dict()
#construct model
from skimage import measure
from scipy.interpolate import griddata
figure(figsize(5,5))
imshow(sumimg)
for mkey in mimgs.keys():
    img = mimgs[mkey]
    contours = [measure.find_contours(img,1)[0]]
    for n, contour in enumerate(contours):
        clen = len(contour[:,0])
        x = griddata(np.arange(clen),contour[:,0],np.linspace(clen,20))
        x = x[~isnan(x)]
        x = hstack((x,x[0]))
        y = griddata(np.arange(clen),contour[:,1],np.linspace(clen,20))
        y = y[~isnan(y)]
        y = hstack((y,y[0]))
        model_data[mkey] = vstack((y,x))
        plot(y, x,linewidth=0.5,color = 'w')
gca().set_xbound((0,1024))
gca().set_ybound((0,1024))
display(gcf());close()
In [8]:
class Basis(dict):    
    def __setitem__(self,key,item):
        """overload the __setitem__ method of dict so the transform and inverse
         transform matrices are computed when the basis vectors are changed"""
        try:
            if key in ['a1','a2']:
                dict.__setitem__(self,key,item)
                A = np.vstack((self['a1'],self['a2'])).T
                A_inv = np.linalg.inv(A)
                self['A'] = A
                self['A_inv'] = A_inv
            else:
                dict.__setitem__(self,key,item)
        except KeyError:
            dict.__setitem__(self,key,item)
                
        
class GeometricModel(object):   
    def __init__(self,lines,basis):
        self.lines = lines
        self.basis = basis
        ## put lines in barycentric coords
        self.barycentric = dict()
        for key in self.lines.keys():
            coords = dot(self.basis['A_inv'],(self.lines[key]-self.basis['p'][:,newaxis])) 
            self.barycentric[key] = coords.T
            
    def coords_from_basis(self,basis):
        ret = dict()
        for key in self.barycentric.keys():
            coords = np.dot(basis['A'],(self.barycentric[key]).T)+basis['p'][:,newaxis]
            ret[key] = coords
        return(ret)

class ModelView(object):
    def __init__(self,model):
        self.model = model
        
    def plot(self,basis,**kwargs):
        lines = self.model.coords_from_basis(basis)
        plot_args = {}
        plot_args['plot_frame'] = kwargs.pop('plot_frame',True)
        plot_args['frame_head_width'] = kwargs.pop('frame_head_width',20)
        plot_args['contour_color'] = kwargs.pop('contour_color','w')
        
        kwargs['color'] = plot_args['contour_color']
        for line in lines.values():
            plot(line[0,:],line[1,:], **kwargs)
        if plot_args['plot_frame']:
            p = basis['p']
            a1 = basis['a1']
            a2 = basis['a2']
            kwargs['color'] = 'g'
            kwargs['head_width'] = plot_args['frame_head_width']
            arrow(p[0],p[1],a1[0],a1[1],**kwargs)
            kwargs['color'] = 'b'
            kwargs['head_width'] = plot_args['frame_head_width']
            arrow(p[0],p[1],a2[0],a2[1],**kwargs)

Now I construct a model and plot it in the reference frame of the confocal image: the position of the three bristles that form the reference frame are hard coded here.

In [9]:
figure(figsize = (10,10))
imfile = tifffile.TiffFile('im_overlay.tiff')
sumimg = imfile.asarray()

###add position of large setae 
model_data['e1'] = np.array([[ 170.02788104,  326.71685254],
                             [ 380.98203222,  919.92627014]])
model_data['e2'] = array([[ 172.83333333,  332.83333333],
                          [ 551.5       ,  164.83333333]])
e1 = model_data['e1']
e2 = model_data['e2']
muscles = dict()

for key in model_data.keys():
    if not(key in ['e1','e2']):
        muscles[key] = model_data[key]
        
confocal_basis = Basis()
confocal_basis['a1'] = e2[1]-e2[0]
confocal_basis['a2'] = e1[1]-e2[0]
confocal_basis['p'] = e2[0]

thorax = GeometricModel(muscles,confocal_basis)
thorax_view = ModelView(thorax)

imshow(sumimg)
import copy
thorax_view.plot(thorax.basis)
gca().axis('tight')
display(gcf());close()

Using the imageAnalysis.py program I manually fit the reference vectors to the gcamp imaging stacks. With the extracted point corespondences, the images from the gcamp movies can then be warped back into the same frame as the confocal stack. Warp and save will perform this on a given fly. I then plot the max projection images from a squadron of eight flies.

In [10]:
# define boolean masks for each muscle
b1_mask = b1_img >0
b2_mask = b2_img >0
b3_mask = b3_img >0

i1_mask = i1_img >0
i2_mask = i2_img >0

iii1_mask = iii1_img >0
iii3_mask = iii3_img >0
iii24_mask = iii24_img >0

hg1_mask = hg1_img >0
hg2_mask = hg2_img >0
hg3_mask = hg3_img >0
hg4_mask = hg4_img >0

tpv_mask = tpv_img >0
tpd_mask = tpd_img >0
pr_mask = pr_img > 0
#Masks for each line
mU56_X_J133 = b1_mask | b2_mask | b3_mask \
            | i1_mask | i2_mask \
            | iii1_mask | iii24_mask |iii3_mask \
            | hg1_mask | hg2_mask | hg3_mask | hg4_mask \
            | tpd_mask | tpv_mask \
            | pr_mask
            
mU82_X_J18 = b1_mask | b2_mask | b3_mask \
            | i1_mask | i2_mask \
            | iii1_mask  |iii3_mask \
            | hg2_mask | hg3_mask | hg4_mask \
            
mU82_X_J100 = iii24_mask | hg2_mask

mU82_X_J22 = b1_mask | b2_mask | b3_mask \
            |iii3_mask \
            | hg1_mask | hg2_mask | hg3_mask | hg4_mask \
            
mU82_X_J23 = b1_mask | b2_mask | b3_mask \
            | iii1_mask |iii3_mask \
            | hg2_mask | hg4_mask \
        
mU82_X_J121 = tpd_mask | tpv_mask
In [11]:
def blur_by_depth(stack):
    from scipy.ndimage.filters import gaussian_filter
    wf = lambda x: (1.0/(x+30))+1.0
    return np.array([gaussian_filter(stack[i],i/2.0 + 3.0)*wf(i) for i in range(shape(stack)[0])])

def model_muscle(stack):
    return sum(blur_by_depth((stack>0).astype(float)),axis = 0)

b1_model = model_muscle(b1_stack)
b2_model = model_muscle(b2_stack)
b3_model = model_muscle(b3_stack)

i1_model = model_muscle(i1_stack)
i2_model = model_muscle(i2_stack)

iii1_model = model_muscle(iii1_stack)
iii3_model = model_muscle(iii3_stack)
iii24_model = model_muscle(iii24_stack)

hg1_model = model_muscle(hg1_stack)
hg2_model = model_muscle(hg2_stack)
hg3_model = model_muscle(hg3_stack)
hg4_model = model_muscle(hg4_stack)

tpd_model = model_muscle(tpd_stack)
tpv_model = model_muscle(tpv_stack)
pr_model = model_muscle(pr_stack)

\(X\) is an matrix of \(m\) pixels by \(n\) frames constructed from the signals from \(k\) muscles so that: \(X = WB\) where \(W\) is a \(m\) by \(k\) mixing matrix and \(B\) is a \(k\) by \(n\) matrix of muscle fluorescence over time. We find the best fit to \(B\) as \(W^\dagger X = B\).

In [12]:
import cv2
from scipy.ndimage.morphology import binary_dilation
bk = binary_dilation(mU56_X_J133,iterations = 60).astype(float)
from scipy.ndimage.filters import gaussian_filter
bk = gaussian_filter(bk,50)
mlist = [b1_model,b2_model,b3_model,
         i1_model,i2_model,
         iii1_model,iii3_model,iii24_model,
         hg1_model,hg2_model,hg3_model,hg4_model,
         tpd_model,tpv_model,pr_model]#,bk]

model_mtrx = [cv2.pyrDown(cv2.pyrDown(ar)) for ar in mlist]
W = np.array([m.ravel() for m in model_mtrx]).T
In [13]:
W_ = W#= W/np.max(W,axis = 0)
imshow(sum(W_.reshape(256,256,-1),axis = 2),cmap = cm.gray)
display(gcf());close()
In [14]:
#Groups for screen
import flylib as flb
import db_access as dba
import cv2
import cPickle
fly_db = dba.get_db()

U82_X_J121 = flb.Squadron(fly_db,[271,272,273,274,275,276])
U82_X_J18 = flb.Squadron(fly_db,[265,267,268,269,270])##266
U82_X_J100 = flb.Squadron(fly_db,[259,260,261,262,263,264])
U82_X_J22 = flb.Squadron(fly_db,[255,256,257,258])
U82_X_J23 = flb.Squadron(fly_db,[278,279,280,281,282]) ##277
U56_X_J133 = flb.Squadron(fly_db,[244,243,242,241,240,239,238,237])

swarm_dict = {'U56_X_J133':U56_X_J133, 
              'U82_X_J18':U82_X_J18, 
              'U82_X_J100':U82_X_J100, 
              'U82_X_J22':U82_X_J22,  
              'U82_X_J23':U82_X_J23,
              'U82_X_J121':U82_X_J121}
import h5py
testfly = swarm_dict['U56_X_J133'].flies[4]
f = h5py.File(testfly.fly_path + "warped_imgs.hdf5", "r")
X_warped = np.array(f['warped_imgs'])
X = X_warped.T.astype(float)
print shape(X)
print shape(W_)
print numpy.linalg.matrix_rank(W_)
W_pin = numpy.linalg.pinv(W_)
B = dot(W_pin,X)
(65536, 35000)
(65536, 15)
15

In [17]:
lns = plot(B.T)
gca().set_xbound(500,2500)
display(gcf());close()
In [20]:
mlist = [b1_model,b2_model,b3_model,
         i1_model,i2_model,
         iii1_model,iii3_model,iii24_model,
         hg1_model,hg2_model,hg3_model,hg4_model,
         tpd_model,tpv_model,pr_model,bk]

model_mtrx = [cv2.pyrDown(cv2.pyrDown(ar)) for ar in mlist]
W = np.array([m.ravel() for m in model_mtrx]).T
W_ = W
W_pin = numpy.linalg.pinv(W_)
B = dot(W_pin,X)
In [25]:
lns = plot(B.T)
gca().set_xbound(500,2500)
gca().set_ybound(-0.5,1.5)
display(gcf());close()
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In [16]:
X_est = dot(W_,B[:,:5000])
#res = X[:,:5000] - X_est
savetiff = hstack((transpose(X_warped[:5000,:].reshape(5000,256,256),(0,2,1))*5,
                   transpose(X_est.T.reshape(5000,256,256),(0,2,1))*5))
tifffile.imsave('model_fit4.tiff',transpose(savetiff.astype(uint8),(0,2,1)))
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In [16]:
savetiff = hstack((transpose(X_warped[:5000,:].reshape(5000,256,256),(0,2,1))*5,
                   transpose(X_est.T.reshape(5000,256,256),(0,2,1))*5))
#savetiff += np.min(savetiff)
#savetiff[savetiff > 256] = 256
tifffile.imsave('model_fit5.tiff',transpose(savetiff.astype(uint8),(0,2,1)))
In []:
 
In []:
 
In [61]:
hist(X_est.T[:5000,:].ravel())
Out[61]:
(array([  2.69167069e+08,   3.56783640e+07,   9.28120300e+06,
          4.33970900e+06,   3.84594400e+06,   3.13025700e+06,
          1.86604100e+06,   3.60070000e+05,   8.71800000e+03,
          2.62500000e+03]),
 array([ -40.60474958,    5.34100217,   51.28675393,   97.23250569,
         143.17825745,  189.12400921,  235.06976097,  281.01551273,
         326.96126449,  372.90701625,  418.85276801]),
 <a list of 10 Patch objects>)
In []:
max(savetiff.ravel())
In [86]:
for i in range(16):
    subplot(16,1,i+1,sharex = gca(),sharey = gca())
    plot(B[i,:])
In [51]:
res = X[:,:5000] - X_est
In [56]:
imshow(res[:,0].reshape(256,256),cmap = cm.gray)
Out[56]:
<matplotlib.image.AxesImage at 0x18c8c42d0>
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In [10]:
import db_access as dba
import flylib as flb
import cv2
import cPickle
output_shape = (1024,1024)
figure(figsize = (8,16))
fly_db = dba.get_db()
swarm = flb.Squadron(fly_db,[244,243,242,241,240,239,238,237])
for i,fly in enumerate(swarm.flies):
    #fly = swarm.flies[0]
    tiffdata = np.array(fly.fly_record['experiments'].values()[0]['tiff_data']['images'])
    pkname = fly.fly_path + '/basis_fits.cpkl'
    #print 'here'
    f = open(pkname)
    img_data = cPickle.load(f)
    f.close()
    test_basis = Basis()
    for key in ['A','p','a1','a2']:
        test_basis[key] = img_data[key]

    src_p0 = test_basis['a1'] + test_basis['p']
    src_p1 = test_basis['p']
    src_p2 = test_basis['a2'] + test_basis['p']
    dst_p0 = confocal_basis['a1'] + confocal_basis['p']
    dst_p1 = confocal_basis['p']
    dst_p2 = confocal_basis['a2'] + confocal_basis['p']
    maximg = np.max(tiffdata,axis=0)
    col = mod(i,2)
    row = i/2
    subplot2grid((4,2),(row,col))
    A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
    imshow(cv2.warpAffine(maximg.T,A,output_shape),cmap = cm.gray)
    thorax_view.plot(thorax.basis,alpha = 0.5)
    gca().set_ybound(0,1024)
    gca().set_xbound(0,1024)
    gca().set_xticks([])
    gca().set_yticks([])
display(gcf());close()

Now I start with a preliminary stab at group analysis.

In [11]:
#Groups for screen
U82_X_J121 = flb.Squadron(fly_db,[271,272,273,274,275,276])
U82_X_J18 = flb.Squadron(fly_db,[265,267,268,269,270])##266
U82_X_J100 = flb.Squadron(fly_db,[259,260,261,262,263,264])
U82_X_J22 = flb.Squadron(fly_db,[255,256,257,258])
U82_X_J23 = flb.Squadron(fly_db,[278,279,280,281,282]) ##277
U56_X_J133 = flb.Squadron(fly_db,[244,243,242,241,240,239,238,237])

# define boolean masks for each muscle
b1_mask = b1_img >0
b2_mask = b2_img >0
b3_mask = b3_img >0

i1_mask = i1_img >0
i2_mask = i2_img >0

iii1_mask = iii1_img >0
iii3_mask = iii3_img >0
iii24_mask = iii24_img >0

hg1_mask = hg1_img >0
hg2_mask = hg2_img >0
hg3_mask = hg3_img >0
hg4_mask = hg4_img >0

tpv_mask = tpv_img >0
tpd_mask = tpd_img >0

#Masks for each line
mU56_X_J133 = b1_mask | b2_mask | b3_mask \
            | i1_mask | i2_mask \
            | iii1_mask | iii24_mask |iii3_mask \
            | hg1_mask | hg2_mask | hg3_mask | hg4_mask \
            | tpd_mask | tpv_mask
            
mU82_X_J18 = b1_mask | b2_mask | b3_mask \
            | i1_mask | i2_mask \
            | iii1_mask  |iii3_mask \
            | hg2_mask | hg3_mask | hg4_mask \
            
mU82_X_J100 = iii24_mask | hg2_mask

mU82_X_J22 = b1_mask | b2_mask | b3_mask \
            |iii3_mask \
            | hg1_mask | hg2_mask | hg3_mask | hg4_mask \
            
mU82_X_J23 = b1_mask | b2_mask | b3_mask \
            | iii1_mask |iii3_mask \
            | hg2_mask | hg4_mask \
        
mU82_X_J121 = tpd_mask | tpv_mask
        
#for plotting
little_frame = copy.copy(thorax.basis)
little_frame['p'] = thorax.basis['p']/4
little_frame['a1'] = thorax.basis['a1']/4
little_frame['a2'] = thorax.basis['a2']/4
In [12]:
i1_select = i1_mask & ~(b2_mask | b3_mask | b1_mask)
i2_select = i2_mask & ~(iii3_mask | hg4_mask | hg1_mask)

b1_select = b1_mask & ~ (i1_mask | b2_mask | b3_mask)
b2_select = b2_mask & ~ (i1_mask | b1_mask | b3_mask | tpd_mask)
b3_select = b3_mask & ~ (i1_mask | b2_mask | b1_mask | tpd_mask)

iii1_select = iii1_mask & ~ (b2_mask | tpd_mask | tpv_mask | iii24_mask)
iii3_select = iii3_mask & ~ (hg1_mask | hg4_mask | hg3_mask | i2_mask)
iii23_select = iii24_mask & ~ (hg1_mask | hg4_mask | hg3_mask | i2_mask | iii3_mask)

hg1_select = hg1_mask & ~ (hg4_mask | hg3_mask | iii3_mask | iii24_mask)
hg2_select = hg2_mask & ~ (hg4_mask)
hg3_select = hg3_mask & ~ (hg2_mask)
hg4_select = hg4_mask & ~ (iii3_mask | iii24_mask)

tp_select = (tpd_mask | tpv_mask) & ~ (i1_mask | b3_mask | iii1_mask | i2_mask)
figure(figsize(12,4))
for i,mask in enumerate([b1_select,b2_select,b3_select,
                         i1_select,i2_select,
                         iii1_select, iii3_select,iii23_select,
                         hg1_select,hg2_select,hg3_select,hg4_select,
                         tp_select]):
    subplot(2,7,i+1)
    imshow(mask)
    thorax_view.plot(thorax.basis,lw = 1.0,plot_frame = False,contour_color = 'white',alpha = 0.5)
    gca().set_xticklabels([]);gca().set_yticklabels([]);
    gca().set_xbound((0,1024));gca().set_ybound((0,1024))
display(gcf());close()
In [19]:
figure()
imshow((hg4_model)+b1_model,cmap=cm.gray)
Out[19]:
<matplotlib.image.AxesImage at 0x16c15c50>
In [21]:
imshow(sum(model_mtrx,axis = 0),cmap = cm.gray)
Out[21]:
<matplotlib.image.AxesImage at 0x8b74890>
In [24]:
plot(B.T)
Out[24]:
[<matplotlib.lines.Line2D at 0x8919310>,
 <matplotlib.lines.Line2D at 0x8919590>,
 <matplotlib.lines.Line2D at 0x89197d0>,
 <matplotlib.lines.Line2D at 0x8919990>,
 <matplotlib.lines.Line2D at 0x8919b50>,
 <matplotlib.lines.Line2D at 0x8919d10>,
 <matplotlib.lines.Line2D at 0x8919ed0>,
 <matplotlib.lines.Line2D at 0x890bb50>,
 <matplotlib.lines.Line2D at 0x891d290>,
 <matplotlib.lines.Line2D at 0x891d450>,
 <matplotlib.lines.Line2D at 0x891d610>,
 <matplotlib.lines.Line2D at 0x891d7d0>,
 <matplotlib.lines.Line2D at 0x891d990>,
 <matplotlib.lines.Line2D at 0x891db50>,
 <matplotlib.lines.Line2D at 0x891d0d0>]
In [98]:
imshow(X_est.T.reshape(-1,256,256)[0],cmap = cm.gray)
figure()
imshow(sum(res.T.reshape(-1,256,256),axis = 0),cmap = cm.gray)
Out[98]:
<matplotlib.figure.Figure at 0x737067d0>
In []:
 
In [14]:
figure(figsize(8,10))
fly = U56_X_J133.flies[1]
expmnt = fly.experiments.values()[0]
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])
rwing = rad2deg(np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])/5.0)
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])
imgs = tifffile.TiffFile(fly.fly_path + 'warped_imgs.tiff').asarray()

ax = subplot2grid((14,6),(0,0),colspan = 5)
plot(times[:30000],ypos[:30000],rasterized=True)
locator_params(axis = 'y', nbins = 3)

subplot2grid((14,6),(1,0),colspan = 5,sharex = ax)
plot(times[:30000],rwing[:30000],rasterized=True)
locator_params(axis = 'y', nbins = 3)
gca().set_ybound(45,90)

for i,mask in enumerate([b1_select,b2_select,b3_select,
                         i1_select,i2_select,
                         iii1_select,iii3_select,iii23_select,
                         hg1_select,hg2_select,hg3_select,hg4_select]):
    subplot2grid((14,6),(2+i,0),colspan = 5,sharex = ax)
    mmask = cv2.pyrDown(cv2.pyrDown((mask).astype(float)))>0
    plot(times[:30000],mean(imgs[:,mmask],axis = 1),rasterized=True)
    if i<11:
        gca().set_xticklabels([])
    locator_params(axis = 'y', nbins = 3)
    
    subplot2grid((14,6),(2+i,5))
    imshow(mask)
    thorax_view.plot(thorax.basis,lw = 1.0,plot_frame = False,contour_color = 'white',alpha = 0.5)
    gca().set_xticklabels([]);gca().set_yticklabels([]);
    gca().set_xbound((0,1024));gca().set_ybound((0,1024))

subplots_adjust(hspace = 0.1)
subplots_adjust(bottom = 0.05)
subplots_adjust(top = 0.99)
subplots_adjust(right = 0.99)
subplots_adjust(left = 0.05)
subplots_adjust(wspace = 0.00)
ax.set_xbound(10,140)
display(gcf());close()
In [90]:
 
In []:
 
In []:
 
In [97]:
from scipy.optimize import nnls
blist = list()
for i in range(1000):
    blist.append(nnls(W,X[:,i]))
In [98]:
B_nn = array([x[0] for x in blist])
In [99]:
X_nn_est = dot(W,B_nn.T)

res_nn = X[:,:1000] - X_nn_est
res = X[:,:1000] - X_est
In [100]:
imshow(X_nn_est.T.reshape(-1,256,256)[0],cmap = cm.gray)
figure()
imshow(sum(res_nn.T.reshape(-1,256,256),axis = 0),cmap = cm.gray)
figure()
imshow(X_est.T.reshape(-1,256,256)[0],cmap = cm.gray)
figure()
imshow(res.T.reshape(-1,256,256)[0],cmap = cm.gray)
Out[100]:
<matplotlib.image.AxesImage at 0x7f864ea6e410>
In [104]:
plot(B.T)
Out[104]:
[<matplotlib.lines.Line2D at 0x7f858f2a7d90>,
 <matplotlib.lines.Line2D at 0x7f858ed2a050>,
 <matplotlib.lines.Line2D at 0x7f858ed2a290>,
 <matplotlib.lines.Line2D at 0x7f858ed2a450>,
 <matplotlib.lines.Line2D at 0x7f858ed2a610>,
 <matplotlib.lines.Line2D at 0x7f858ed2a7d0>,
 <matplotlib.lines.Line2D at 0x7f858ed2a990>,
 <matplotlib.lines.Line2D at 0x7f858f296f10>,
 <matplotlib.lines.Line2D at 0x7f858ed2ad10>,
 <matplotlib.lines.Line2D at 0x7f858ed2aed0>,
 <matplotlib.lines.Line2D at 0x7f858ed2c0d0>,
 <matplotlib.lines.Line2D at 0x7f858ed2c290>,
 <matplotlib.lines.Line2D at 0x7f858ed2c450>,
 <matplotlib.lines.Line2D at 0x7f858ed2c610>]
In [103]:
close()
In []:
 
In []:
 
In []:
 
In []:
 
In [94]:
#this block will load the metadata and warp the images/movies into a common reference frame.
swarm_dict = {'U56_X_J133':U56_X_J133, 
              'U82_X_J18':U82_X_J18, 
              'U82_X_J100':U82_X_J100, 
              'U82_X_J22':U82_X_J22,  
              'U82_X_J23':U82_X_J23,
              'U82_X_J121':U82_X_J121}

def warp_fly_movie(fly):
    fname = fly.fly_path + 'meta_data.cpkl'
    f = open(fname,'rb');data = cPickle.load(f);f.close()
    pkname = fly.fly_path + '/basis_fits.cpkl'
    f = open(pkname);img_data = cPickle.load(f);f.close()
    test_basis = Basis()
    #print data.keys()
    for key in ['A','p','a1','a2']:
        test_basis[key] = img_data[key]
    src_p0 = test_basis['a1'] + test_basis['p']
    src_p1 = test_basis['p']
    src_p2 = test_basis['a2'] + test_basis['p']
    dst_p0 = confocal_basis['a1'] + confocal_basis['p']
    dst_p1 = confocal_basis['p']
    dst_p2 = confocal_basis['a2'] + confocal_basis['p']
    A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
    #compose back to Ca image size
    A = vstack((A,[0,0,1]))
    S = np.array([[0.25,0,0],[0,0.25,0],[0,0,1]])
    Ap = dot(S,A)[:2,:3]
    testmovie = np.array(fly.fly_record['experiments'].values()[0]['tiff_data']['images'])
    X_warped = [cv2.warpAffine(testmovie[frame].T,Ap,(256,256)).ravel() for frame in range(shape(testmovie)[0])]
    X_warped = np.array(X_warped)
    import h5py
    f = h5py.File(fly.fly_path + "warped_imgs.hdf5", "a")
    f['warped_imgs'] = X_warped
    f.flush()
    f.close()
In [95]:
for swarm in swarm_dict.values():
    for fly in swarm.flies:
#        os.remove(fly.fly_path + "warped_imgs.hdf5")
        print fly.fly_num
        warp_fly_movie(fly)
255
256
257
258
278
279
280
281
282
271
272
273
274
275
276
259
260
261
262
263
264
265
267
268
269
270
244
243
242
241
240
239
238
237

In [87]:
plot(B[:,:5000].T)
Out[87]:
[<matplotlib.lines.Line2D at 0x17e513dd0>,
 <matplotlib.lines.Line2D at 0x17e518090>,
 <matplotlib.lines.Line2D at 0x17e5182d0>,
 <matplotlib.lines.Line2D at 0x17e518490>,
 <matplotlib.lines.Line2D at 0x17e518650>,
 <matplotlib.lines.Line2D at 0x17e518810>,
 <matplotlib.lines.Line2D at 0x17e5189d0>,
 <matplotlib.lines.Line2D at 0x1bf430150>,
 <matplotlib.lines.Line2D at 0x17e518d50>,
 <matplotlib.lines.Line2D at 0x17e518f10>,
 <matplotlib.lines.Line2D at 0x17e519110>,
 <matplotlib.lines.Line2D at 0x17e5192d0>,
 <matplotlib.lines.Line2D at 0x17e519490>,
 <matplotlib.lines.Line2D at 0x17e519650>]
In [88]:
X_est = dot(W,B[:,:5000])
In [58]:
print mean(B[:,:5000],axis = 1)
[ 1.59878007  0.49757645  1.77607475  0.43389707  1.21200045  0.34361473
  1.63750865  0.58052898  0.11960622  1.1447411   0.62986823 -0.65435049
  0.55309889  0.62362939]

In [60]:
from scipy.optimize import nnls
In [133]:
 
Out[133]:
(5000, 512, 256)
In [56]:
imshow(Out[51][0])
Out[56]:
<matplotlib.image.AxesImage at 0x3f701dd0>
In [15]:
#this block will load the metadata and warp the images/movies into a common reference frame.
swarm_dict = {'U56_X_J133':U56_X_J133, 
              'U82_X_J18':U82_X_J18, 
              'U82_X_J100':U82_X_J100, 
              'U82_X_J22':U82_X_J22,  
              'U82_X_J23':U82_X_J23,
              'U82_X_J121':U82_X_J121}



img_dict = dict()
output_shape = (256,256)
import cPickle
for swarm_key in swarm_dict.keys():
    swarm = swarm_dict[swarm_key]
    
    corr_list = list()
    dot_list = list()
    max_list = list()
    wing_list = list()
    norm_list = list()
    mean_list = list()

    movie_list = list()
    rwing_list = list()
    ypos_list = list()

    for fly in swarm.flies:
        try:
            #print(fly.fly_num)
            fname = fly.fly_path + 'meta_data.cpkl'
            #f = open(fname,'rb');data = cPickle.load(f)['cl_stim_metadata'];f.close()
            f = open(fname,'rb');data = cPickle.load(f);f.close()
            pkname = fly.fly_path + '/basis_fits.cpkl'
            f = open(pkname);img_data = cPickle.load(f);f.close()
            test_basis = Basis()
            #print data.keys()
            for key in ['A','p','a1','a2']:
                test_basis[key] = img_data[key]
            src_p0 = test_basis['a1'] + test_basis['p']
            src_p1 = test_basis['p']
            src_p2 = test_basis['a2'] + test_basis['p']
            dst_p0 = confocal_basis['a1'] + confocal_basis['p']
            dst_p1 = confocal_basis['p']
            dst_p2 = confocal_basis['a2'] + confocal_basis['p']
            A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
            #compose back to Ca image size
            A = vstack((A,[0,0,1]))
            S = np.array([[0.25,0,0],[0,0.25,0],[0,0,1]])
            Ap = dot(S,A)[:2,:3]
            #warp the images
            dot_list.append(cv2.warpAffine(data['sin_stim_metadata']['dot_img'].T,Ap,output_shape))
            corr_list.append(cv2.warpAffine(data['sin_stim_metadata']['corr_img'].T,Ap,output_shape))
            max_list.append(cv2.warpAffine(data['sin_stim_metadata']['max_img'].T,Ap,output_shape))
            mean_list.append(cv2.warpAffine(data['sin_stim_metadata']['mean_img'].T,Ap,output_shape))
            norm_list.append(cv2.warpAffine(data['sin_stim_metadata']['norm_img'].T,Ap,output_shape))
            #now warp the movies
            new_shape = (shape(data['ave_movie'])[0],output_shape[0],output_shape[1])
            warped_ave_movie = zeros(new_shape)
            for frame in arange(shape(warped_ave_movie)[0]):
                warped_ave_movie[frame,:,:] = cv2.warpAffine(data['ave_movie'][frame,:,:].T,Ap,output_shape)
            movie_list.append(warped_ave_movie/mean_list[-1])
            rwing_list.append(data['ave_rwing'])
            ypos_list.append(data['ave_ypos'])
        except Exception as inst:
            print inst
        img_dict[swarm_key] = {'dot_list':dot_list,
                                'corr_list':corr_list,
                                'mean_list':mean_list,
                                'norm_list':norm_list,
                                'max_list':max_list,
                                'movie_list':movie_list,
                                'rwing_list':rwing_list,
                                'ypos_list':ypos_list}
-c:62: RuntimeWarning: invalid value encountered in divide

In [16]:
ndat = np.array(img_dict[swarm_key]['norm_list'][:]).ravel()
ddat = np.array(img_dict[swarm_key]['dot_list'][:]).ravel()
cdat = np.array(img_dict[swarm_key]['corr_list'][:]).ravel()

norm_max = np.percentile(ndat,95)
dot_max = np.percentile(ddat,95)
corr_max = np.percentile(cdat,95)

norm_min = np.percentile(ndat,0)
dot_min = np.percentile(ddat,5)
corr_min = np.percentile(cdat,5)

def summarize_swarm(swarm_key,swarm_mask,axlist):
    axes(axlist[0])
    imshow(swarm_mask)
    thorax_view.plot(thorax.basis,lw = 1.0,plot_frame = False)
    axes(axlist[1])
    imshow(np.mean(img_dict[swarm_key]['norm_list'],axis = 0),cmap = cm.gray,vmin = 0, vmax = norm_max)
    thorax_view.plot(little_frame,lw = 1.0,plot_frame = False,contour_color = 'orange',alpha = 0.5)
    axes(axlist[2])
    imshow(np.mean(img_dict[swarm_key]['dot_list'],axis = 0),cmap = cm.gray,vmin = dot_min, vmax = dot_max)
    thorax_view.plot(little_frame,lw = 1.0,plot_frame = False,contour_color = 'orange',alpha = 0.5)
    axes(axlist[3])
    imshow(np.mean(img_dict[swarm_key]['corr_list'],axis = 0),cmap = cm.gray,vmin = corr_min, vmax = corr_max)
    thorax_view.plot(little_frame,lw = 1.0,plot_frame = False,contour_color = 'orange',alpha = 0.5)

figure(figsize(8,10))
axlist = [subplot(5,4,i) for i in range(1,5)]
summarize_swarm('U56_X_J133',mU56_X_J133,axlist)
axlist[0].set_ylabel("U56 X J133")

#axlist[1].set_title(r"$\sqrt{diag(XX^\top)}$")
#axlist[2].set_title(r"$Xs^\top$")
#axlist[3].set_title(r"$\frac{Xs^\top}{\sqrt{diag(XX^\top)}||s||}$",y= 1.08)

axlist[1].set_title(r"pixel std")
axlist[2].set_title(r"dot image")
axlist[3].set_title(r"correlation image")

axlist = [subplot(5,4,i) for i in range(5,9)]
summarize_swarm('U82_X_J18',mU82_X_J18,axlist)
axlist[0].set_ylabel("U82 X J18")

axlist = [subplot(5,4,i) for i in range(9,13)]
summarize_swarm('U82_X_J100', mU82_X_J100,axlist)
axlist[0].set_ylabel("U82 X J100")

axlist = [subplot(5,4,i) for i in range(13,17)]
summarize_swarm('U82_X_J22', mU82_X_J22,axlist)
axlist[0].set_ylabel("U82 X J22")

axlist = [subplot(5,4,i) for i in range(17,21)]
summarize_swarm('U82_X_J23', mU82_X_J23,axlist)
axlist[0].set_ylabel("U82 X J23")
for ax in gcf().axes:
    ax.set_xticklabels([])
    ax.set_yticklabels([])
    ax.set_xticks([]); ax.set_yticks([])
for i in [1,2,3,
          5,6,7,
          9,10,11,
          13,14,15,
          17,18,19]:
    ax = gcf().axes[i]
    ax.set_xbound(15,190);ax.set_ybound(60,236)
for i in [0,4,8,12,16]:
    ax = gcf().axes[i]
    ax.set_xbound(15*4,190*4);ax.set_ybound(60*4,235*4)
#tight_layout()
display(gcf());#close()
In [17]:
swarm_dict = {'U56_X_J133':U56_X_J133, 
              'U82_X_J18':U82_X_J18, 
              'U82_X_J100':U82_X_J100, 
              'U82_X_J22':U82_X_J22,  
              'U82_X_J23':U82_X_J23,
              'U82_X_J121':U82_X_J121}

m_dict = {'U56_X_J133':mU56_X_J133,
         'U82_X_J18':mU82_X_J18,
         'U82_X_J100':mU82_X_J100,
         'U82_X_J22':mU82_X_J22, 
         'U82_X_J23':mU82_X_J23,
         'U82_X_J121':mU82_X_J121}

def plot_group_sig(sig):
    sig_mean = np.mean(sig,axis = 0)
    sig_std = np.std(sig,axis = 0)
    times = np.linspace(0,2,60*2)
    #fill_between(times,sig_mean+sig_std,sig_mean-sig_std,alpha = 0.5)
    plot(times,sig.T,color = 'k',alpha = 0.5)#sig_mean)
    plot(times,sig_mean,lw = 2.0,alpha = 0.8)

figure(figsize = (8,10))
for j,swarm_key in enumerate(swarm_dict.keys()):
    #swarm_key = 'U82_X_J23'
    subplot2grid((13,7),(0,j),colspan = 1)
    imshow(m_dict[swarm_key])
    thorax_view.plot(thorax.basis,lw = 0.5,alpha = 0.5,plot_frame = False)
    gca().set_xbound(0,1024);gca().set_ybound(0,1024)
    
    wing_angle = rad2deg(np.array(img_dict[swarm_key]['rwing_list'])/5)
    subplot2grid((13,7),(1,j),colspan = 1)
    plot_group_sig(wing_angle)
    gca().set_ybound(45,75)

    def plot_roi(mask):
        mmask = cv2.pyrDown(cv2.pyrDown(mask.astype(float)))>0
        mdata = np.mean(np.array(img_dict[swarm_key]['movie_list'])[:,:,mmask],axis=2)
        plot_group_sig(mdata)
        gca().set_ybound(0.9,1.1)

    def imshow_roi(mask):
        mmask = cv2.pyrDown(cv2.pyrDown(mask.astype(float)))>0
        imshow(mmask)
        thorax_view.plot(little_frame,lw = 0.5,plot_frame = False,contour_color = 'white',alpha = 0.5)
        gca().set_xticklabels([]);gca().set_yticklabels([]);
        gca().set_xbound(15,190);gca().set_ybound(60,236)

    for i,mask in enumerate([b1_select,
                             b2_select,
                             b3_select,
                             i1_select,
                             i2_select,
                             iii1_select,
                             iii3_select,
                             hg1_select,
                             hg2_select,
                             hg3_select,
                             hg4_select,
                             tp_select]):
        subplot2grid((14,7),(i+2,j),colspan = 1)
        plot_roi(mask)
        subplot2grid((14,7),(i+2,6))
        imshow_roi(mask)

for ax in gcf().axes:
    ax.set_yticklabels([])
    ax.set_xticklabels([])
    
subplots_adjust(hspace = 0.1)
subplots_adjust(bottom = 0.05)
subplots_adjust(top = 0.99)
subplots_adjust(right = 0.99)
subplots_adjust(left = 0.05)
subplots_adjust(wspace = 0.00)
display(gcf())
In [20]:
print(shape(mask_mtrx))
print(numpy.linalg.matrix_rank(mask_mtrx))
(14, 65536)
14

In [21]:
imshow(sum(mask_mtrx,axis = 0).reshape(256,256))
Out[21]:
<matplotlib.image.AxesImage at 0x3aa5f510>
In [24]:
 
In [42]:
X_mean = np.mean(X_warped.astype(float),axis = 0)
#X_Xbar = X_warped/X_mean
X_norm = np.linalg.norm(X_warped.astype(float),axis = 0)
In [59]:
X_Xnorm = (X_warped.astype(float)-X_mean)[~(xwarped == 0)]/X_norm
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
<ipython-input-59-62a6d978f56b> in <module>()
----> 1 X_Xnorm = (X_warped.astype(float)-X_mean)[~(xwarped == 0)]/X_norm

MemoryError: 
In [58]:
hist(X_Xnorm[~isnan(X_Xnorm)].ravel())
Out[58]:
(array([  2.64702270e+08,   9.08990000e+04,   5.14500000e+03,
          1.38400000e+03,   2.40000000e+02,   4.00000000e+01,
          5.00000000e+00,   1.20000000e+01,   3.00000000e+00,
          2.00000000e+00]),
 array([-0.01414072,  0.05795574,  0.13005221,  0.20214868,  0.27424514,
         0.34634161,  0.41843807,  0.49053454,  0.56263101,  0.63472747,
         0.70682394]),
 <a list of 10 Patch objects>)
In [60]:
del(np)
In [47]:
imshow(X_mean.reshape(256,256))
Out[47]:
<matplotlib.image.AxesImage at 0xb42cf310>
In [213]:
gc.collect()
Out[213]:
4728
In [200]:
plot(X_Xbar[500])
Out[200]:
[<matplotlib.lines.Line2D at 0x7f8c4566c190>]
In [129]:
imshow(fit_movie.T[9].reshape(256,256),cmap = cm.gray)
Out[129]:
<matplotlib.image.AxesImage at 0xb5793490>
In [134]:
tifffile.imsave('fittest2.tiff',fit_movie.T.reshape(10000,256,256).astype(uint8)*10)
In [138]:
tifffile.imsave('fittest_X.tiff',X_warped.reshape(10000,256,256).astype(uint8))
In []:
 
In [85]:
#functions to summarize the sin stimulus protocol
def collect_cycles(ypos,stim_cond,trial_selection_mask,flight_mask):
    cycles = list()
    trial_selector = zeros_like(ypos).astype(bool)
    cycle_selection_mask = ypos>1.5
    #trial_selection_mask = around(stim_cond)==2
    trial_idxs = flb.idx_by_thresh(trial_selection_mask,0.5)
    for trial_idx in trial_idxs:
        if not(np.sum((~flight_mask)[trial_idx])):
            trial_selector[:] = 0
            trial_selector[trial_idx] = 1
            cycle_idxs = flb.idx_by_thresh(trial_selector & cycle_selection_mask)
            cycles.extend([arange(cycle_idxs[i][0],cycle_idxs[i+1][0]) 
                       for i in range(len(cycle_idxs)-1)])
    return cycles

def cycle_average(sig,cycle_idxs,times,flight_mask):
    """function extracts the cycle average of the signal given a list of idxs where the cycles 
    happen, and a mask to exclude non-flight bouts. Resamples so that cycle samples line up.
    returns a dictionary of summary stats"""
    trial_idxs = [np.arange(x[0],x[0]+210) for x in cycle_idxs]# if (len(x)>0 and (sum(~flight_mask[x[0]:x[0]+200])==0))]
    trial_times = [times[idx]-times[idx[0]] for idx in trial_idxs if sum(~flight_mask[idx])==0]
    trial_data = [sig[idx] for idx in trial_idxs if sum(~flight_mask[idx])==0]
    import scipy.interpolate
    fun = scipy.interpolate.griddata
    time_resamp = linspace(0,2.,60*2.)
    data_resamp = list()
    for td,tt in zip(trial_data,trial_times):
        f = scipy.interpolate.interp1d(tt, td, kind='linear', axis=0)
        data_resamp.append(f(time_resamp)) 
    retdict =   {'time_resamp':time_resamp,
                'sig_resamp':data_resamp}
    return retdict

def wing_correlations(sig,imgs,fly = None):
    #calculate some summary images including the correlation image
    mean_img = mean(imgs,axis=0)
    zeroed_img = imgs - mean_img
    max_img = np.max(imgs,axis = 0)
    norm_img = numpy.linalg.norm(zeroed_img,axis = 0)
    norm_sig = numpy.linalg.norm(sig-np.mean(sig))
    dot_img = sum(zeroed_img*(sig[:,newaxis,newaxis]-np.mean(sig)),axis = 0)
    corr_img = dot_img/(norm_img*norm_sig)
    return {'mean_img':mean_img,
        'max_img':max_img,
        'norm_img':norm_img,
        'dot_img':dot_img,
        'corr_img':corr_img}
    
def calc_fly_stats(fly):
    #only did one experiment on these flies
    expmnt = fly.experiments.values()[0]
    nframes = 30000
    #load the movie
    movie = np.array(expmnt.exp_record['tiff_data']['images'][:nframes])
    flight_mask = np.array(expmnt.exp_record['flight_mask'][:nframes])
    #Map some variables
    stim_cond = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['DAC2'])[:nframes]
    lwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph0'])[:nframes]
    rwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])[:nframes]
    lmr = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph2'])[:nframes]
    aux = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph3'])[:nframes]
    xpos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Xpos'])[:nframes]
    ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])[:nframes]
    times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])[:nframes]
    #idxs associated with sin type stimulus
    trial_selection_mask = stim_cond>0
    cycle_idxs = collect_cycles(ypos,stim_cond,trial_selection_mask,flight_mask)
    #concatinate the stimuli idxs
    sin_idxs = hstack(cycle_idxs)
    #select the frames from the imaging stack 
    #that correspond to the sin stimulus
    imgs = movie[sin_idxs,:,:]
    sig = rwing[sin_idxs]
    #correlate with the right wingstroke data
    sin_stim_metadata = wing_correlations(sig,imgs,fly)
    #now do the same thing for the stripe-fixation epoch at the end of the trial
    end_cl_idx = flb.idx_by_thresh(stim_cond<0)[-1]
    cl_selection_mask = np.zeros_like(flight_mask)
    cl_selection_mask[end_cl_idx] = 1
    cl_idxs = squeeze(np.argwhere(cl_selection_mask & flight_mask))
    imgs = movie[cl_idxs,:,:]
    sig = rwing[cl_idxs]
    #again correlate with the right wingstroke data
    cl_stim_metadata = wing_correlations(sig,imgs,fly)
    #now cycle-average the image stream and left and right wing signals
    rdict = cycle_average(movie,cycle_idxs,times,flight_mask)
    ave_movie = np.mean(rdict['sig_resamp'],axis = 0)
    rdict = cycle_average(rwing,cycle_idxs,times,flight_mask)
    ave_rwing = np.mean(rdict['sig_resamp'],axis = 0)
    rdict = cycle_average(lwing,cycle_idxs,times,flight_mask)
    ave_lwing = np.mean(rdict['sig_resamp'],axis = 0)
    rdict = cycle_average(ypos,cycle_idxs,times,flight_mask)
    ave_ypos = np.mean(rdict['sig_resamp'],axis = 0)
    return {'sin_stim_metadata':sin_stim_metadata,
            'cl_stim_metadata':cl_stim_metadata,
            'ave_movie':ave_movie,
            'ave_rwing':ave_rwing,
            'ave_lwing':ave_lwing,
            'ave_ypos':ave_ypos}

#calculates and save the 'metadata' for each fly.
#import cPickle
#for swarm in [U82_X_J121, U82_X_J18, U82_X_J100,U82_X_J22,  U82_X_J23, U56_X_J133]:
#    #for swarm in [U82_X_J18]:
#    for fly in swarm.flies:
#        try:
#            #print(fly.fly_num)
#            rdict = calc_fly_stats(fly)
#            fname = fly.fly_path + 'meta_data.cpkl'
#            f = open(fname,'wb')
#            cPickle.dump(rdict,f)
#            f.close()
#            del(rdict)
#        except Exception as inst:
#            print inst
In [88]:
for swarm_key in swarm_dict.keys():
    print swarm_key
    group_movie = np.mean(img_dict[swarm_key]['movie_list'],axis = 0)
    movie_std = np.std(group_movie[~isnan(group_movie)])
    movie_min = 1-movie_std*8
    movie_max = 1+movie_std*8
    rescaled = (group_movie-movie_min)*255/(movie_std*16)
    rescaled[isnan(rescaled)] = 255/2
    rescaled[rescaled>255] = 255
    rescaled[rescaled<0] = 0
    tifffile.imsave(swarm_key + '.tiff',np.uint8(rescaled))
U82_X_J22
U82_X_J23
U82_X_J121
U82_X_J100
U82_X_J18
U56_X_J133

In [86]:
 
In [83]:
 
In [186]:
 
In [186]:
 
In [186]:
 
In [98]:
for i,p in enumerate([0,30,60,90]):
    subplot(1,4,i+1)
    imshow(group_movie[p],cmap = cm.gray,vmin = movie_min,vmax = movie_max)
In [79]:
 
In [75]:
mov = np.array(img_dict[swarm_key]['movie_list'])
In [103]:
tifffile.imsave('test2.tiff',uint8(mov[1,:,:,:]*50))
In [103]:
imshow(hg3_mask & ~(hg4_mask | i2_mask | hg1_mask))
Out[103]:
<matplotlib.image.AxesImage at 0x1a8a43350>
In [113]:
shape(img_dict[swarm_key]['rwing_list'])
Out[113]:
(5, 120)
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In [19]:
 
In [20]:
 
In [21]:
 
Out[21]:
[<matplotlib.lines.Line2D at 0xc862d10>]
In [32]:
display(gcf())
In [34]:
plot(rwing[25000])
Out[34]:
[<matplotlib.lines.Line2D at 0x3afdced0>]
In [23]:
fm = np.array(expmnt.exp_record['flight_mask'])
In [37]:
argwhere(~fm)[0]
Out[37]:
array([23991])
In [43]:
imshow(expmnt.exp_record['tiff_data']['images'][int(argwhere(~fm))])
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-43-2b13314dfe6b> in <module>()
----> 1 imshow(expmnt.exp_record['tiff_data']['images'][int(argwhere(~fm)[1000])])

IndexError: index 1000 is out of bounds for axis 0 with size 275
In [46]:
np.mean(np.array(expmnt.exp_record['tiff_data']['images'][int(argwhere(~fm))]),axis = 0)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-46-65a02ecb4a58> in <module>()
----> 1 np.mean(np.array(expmnt.exp_record['tiff_data']['images'][int(argwhere(~fm))]),axis = 0)

TypeError: only length-1 arrays can be converted to Python scalars
In [24]:
F = np.mean(np.array(expmnt.exp_record['tiff_data']['images'])[~fm,:,:],axis = 0)
In [25]:
DF = np.array(expmnt.exp_record['tiff_data']['images'])[fm,:,:]
In [26]:
normed = (DF-F)/F
In [34]:
normed[normed<0] = 0
In [37]:
tifffile.imsave('nstream.tiff', uint8(normed[:500]*50))
In [39]:
plot(rwing[fm])
Out[39]:
[<matplotlib.lines.Line2D at 0x360a75d0>]
In [33]:
hist(uint8(normed[:500].ravel()*255))
Out[33]:
(array([ 1530393.,   767034.,   537929.,   436250.,   491659.,   656907.,
         1175996.,  2047392.,  2007581.,  1658859.]),
 array([   0. ,   25.5,   51. ,   76.5,  102. ,  127.5,  153. ,  178.5,
         204. ,  229.5,  255. ]),
 <a list of 10 Patch objects>)
In [165]:
swarm = U82_X_J18
fly = swarm.flies[0]
In [169]:
fname = fly.fly_path + 'meta_data.cpkl'
f = open(fname,'rb');data = cPickle.load(f);f.close()
In [170]:
data.keys()
Out[170]:
['cl_stim_metadata',
 'ave_movie',
 'sin_stim_metadata',
 'ave_rwing',
 'ave_lwing']
In [172]:
pkname = fly.fly_path + '/basis_fits.cpkl'
f = open(pkname);img_data = cPickle.load(f);f.close()
test_basis = Basis()
for key in ['A','p','a1','a2']:
    test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
#compose back to Ca image size
A = vstack((A,[0,0,1]))
S = np.array([[0.25,0,0],[0,0.25,0],[0,0,1]])
Ap = dot(S,A)[:2,:3]
In [174]:
imshow(data['ave_movie'][0,:,:])
Out[174]:
<matplotlib.image.AxesImage at 0x35721950>
In [182]:
 
In [186]:
 
In [190]:
imshow(warped_ave_movie[50,:,:])
Out[190]:
<matplotlib.image.AxesImage at 0x307f1210>
In [163]:
 
In [163]:
 
In [163]:
 
In [163]:
 
In [164]:
imshow(f(np.linspace(0,2,5000)))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-164-fcfd7d3058f6> in <module>()
----> 1 imshow(f(np.linspace(0,2,5000)))

TypeError: 'file' object is not callable
In [321]:
tt
Out[321]:
array([ 0.      ,  0.014105,  0.02821 ,  0.04231 ,  0.05641 ,  0.07051 ,
        0.08461 ,  0.09871 ,  0.112815,  0.12692 ,  0.14102 ,  0.15512 ,
        0.16922 ,  0.18332 ,  0.197425,  0.21153 ,  0.22563 ,  0.23973 ,
        0.25383 ,  0.26793 ,  0.282035,  0.296135,  0.31024 ,  0.32434 ,
        0.33844 ,  0.35254 ,  0.36664 ,  0.380745,  0.39485 ,  0.40895 ,
        0.42305 ,  0.43715 ,  0.45125 ,  0.465355,  0.47946 ,  0.49356 ,
        0.50766 ,  0.52176 ,  0.53586 ,  0.54996 ,  0.564065,  0.57817 ,
        0.59227 ,  0.60637 ,  0.62047 ,  0.63457 ,  0.648675,  0.66278 ,
        0.67688 ,  0.69098 ,  0.70508 ,  0.71918 ,  0.73328 ,  0.747385,
        0.76149 ,  0.77559 ,  0.78969 ,  0.80379 ,  0.81789 ,  0.831995,
        0.8461  ,  0.8602  ,  0.8743  ,  0.8884  ,  0.9025  ,  0.916605,
        0.93071 ,  0.94481 ,  0.95891 ,  0.97301 ,  0.98711 ,  1.00121 ,
        1.015315,  1.02942 ,  1.04352 ,  1.05762 ,  1.07172 ,  1.08582 ,
        1.099925,  1.11403 ,  1.12813 ,  1.14223 ,  1.15633 ,  1.17043 ,
        1.18453 ,  1.198635,  1.21274 ,  1.22684 ,  1.24094 ,  1.25504 ,
        1.26914 ,  1.283245,  1.29735 ,  1.31145 ,  1.32555 ,  1.33965 ,
        1.35375 ,  1.367855,  1.381955,  1.39723 ,  1.41133 ,  1.42543 ,
        1.43953 ,  1.453635,  1.46774 ,  1.48184 ,  1.49594 ,  1.51004 ,
        1.52414 ,  1.538245,  1.55235 ,  1.56645 ,  1.58055 ,  1.59465 ,
        1.60875 ,  1.62285 ,  1.636955,  1.65106 ,  1.66516 ,  1.67926 ,
        1.69336 ,  1.70746 ,  1.721565,  1.73567 ,  1.74977 ,  1.76387 ,
        1.77797 ,  1.79207 ,  1.806175,  1.82028 ,  1.83438 ,  1.84848 ,
        1.86258 ,  1.87668 ,  1.89078 ,  1.904885,  1.91899 ,  1.93309 ,
        1.94719 ,  1.96129 ,  1.97539 ,  1.989495,  2.0036  ])
In [123]:
 
In [123]:
 
In [123]:
 
In [196]:
 
In [197]:
 
Out[197]:
0.74133576559200676
In [198]:
 
In [199]:
 
Out[199]:
-0.72825339992174254
In [201]:
 
Out[201]:
2594.2955387660832
In [123]:
 
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 

To perform ICA I take the first fly in the squadron and warp the movie into a scaled-down version of the confocal frame. The projection of the segemented muscles from the confocal stacks can then be used as pixel masks on the transformed image stack. An alternative would be to warp the masks into the Ca++ imaging frame; the approach I am taking now is conceptually simpler though no doubt computationally much slower. In the analysis these masks are used in two ways: 1) reduce the size of the problem 2) provide spatial information

In []:
 

First I performed ICA on a mask defined as the union of all the muscles. The 12 columns of the mixing matrx are reshaped into annatomical space and shown as a heat map on the left. The corresponding component is shown to the right. Time is in sec.

In [77]:
ica_filters, ica_sig = st_ica(mask_signal, ncomp = 7,  mu = 0.8)
In [89]:
for i in range(shape(ica_filters)[0]):
    tmpimg = zeros((256,256))
    tmpimg[mask] = ica_filters[i,:]+abs(np.min(ica_filters[i,:]))
    figure()
    little_frame = copy.copy(thorax.basis)
    little_frame['p'] = thorax.basis['p']/4
    little_frame['a1'] = thorax.basis['a1']/4
    little_frame['a2'] = thorax.basis['a2']/4
    thorax_view.plot(little_frame,lw = 0.3)
    imshow(tmpimg)
In [16]:
fly =U_X_J.flies[0]
#load the warped movie
warped_movie = tifffile.TiffFile(fly.fly_path + 'warped_imgs.tiff').asarray()
#only did one experiment on these flies
expmnt = fly.experiments.values()[0]
#Map some variables
stim_cond = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['DAC2'])
lwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph0'])
rwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])
lmr = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph2'])
aux = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph3'])
xpos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Xpos'])
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])

flight_mask = aux < 0.2
trial_selection_mask = stim_cond>0
cycle_idxs = collect_cycles(ypos,stim_cond,trial_selection_mask)
In [21]:
nc = 7
posterior_muscle_mask = i2_mask | iii3_mask \
                        | iii24_mask | hg1_mask | hg2_mask \
                        | hg3_mask | hg4_mask 

ica_data = warped_movie[~flight_mask[:5000],:,:]

S_,A_,C_ = runICA(ica_data,posterior_muscle_mask,n_components = nc)

plot_ICA(warped_movie,
        posterior_muscle_mask,
        flight_mask,
        times,
        rwing,
        cycle_idxs,
        S_,A_,C_)

plot_comp_img(A_,posterior_muscle_mask)
/home/flyranch/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/sklearn/decomposition/fastica_.py:110: UserWarning: FastICA did not converge. You might want to increase the number of iterations.
  ' to increase the number of iterations.')

In [15]:
def runICA(warped_movie,mask,n_components = 2):
    """get the ICA signals,mixing matrix and unmixing matrix given movie,
    mask and expected number of components"""
    mask = cv2.pyrDown(cv2.pyrDown(np.float64(mask))) >0
    mask_signal = warped_movie[:,mask]
    mask_signal = np.float64(mask_signal)
    from sklearn.decomposition import FastICA, PCA
    ica = FastICA(n_components=n_components,max_iter=1000)
    S  = mask_signal/mask_signal.std(axis=0)
    S_ = ica.fit_transform(mask_signal)  # Reconstruct signals
    A_ = ica.mixing_  # Get estimated mixing matrix
    C_ = ica.components_
    return S_,A_,C_

def plot_comp_img(A_,mask):
    """just plot a grid of the component images"""
    figure(figsize = (8.5,12))
    mask = cv2.pyrDown(cv2.pyrDown(mask.astype(float))) > 0
    n_components = shape(A_)[1]
    n_rows = int(ceil(nc/3.0))
    for i in range(n_components):
        col = mod(i,3)
        row = i/3
        subplot2grid((n_rows,3),(row,col)) 
        tst_img = np.zeros((256,256))
        tst_img[mask] = A_[:,i] - np.min(A_)
        imshow(tst_img,vmin= 0)#,vmax = 100)
        little_frame = copy.copy(thorax.basis)
        little_frame['p'] = thorax.basis['p']/4
        little_frame['a1'] = thorax.basis['a1']/4
        little_frame['a2'] = thorax.basis['a2']/4
        thorax_view.plot(little_frame,lw = 0.3)
        gca().set_xticks([]);gca().set_yticks([])
        gca().set_xbound((0,256));gca().set_ybound((0,256))
    
def plot_ICA(warped_movie,
             mask,
             flight_mask,
             times,
             rwing,
             cycle_idxs,
             S_,A_,C_):
    """plot a summary of the ICA data from the entire experiment using the components
    extracted from a small subset"""
    n_components = shape(A_)[1]
    n_samples = sum(~flight_mask[:30000])
    mask = cv2.pyrDown(cv2.pyrDown(mask.astype(float))) > 0
    dat = warped_movie[~flight_mask[:30000],:,:][:,mask]
    sigs = np.dot(C_,dat.T)
    figure(figsize = (8.5,12))
    rwing = rad2deg(rwing/5.0)
    #plot the rwing amp cycle average
    subplot2grid((n_components+2,5),(0,4))
    cd = cycle_average(rwing,cycle_idxs,times)
    plot(cd['time_resamp'],cd['sig_mean'])
    #plot the rwing example trace
    subplot2grid((n_components+2,5),(0,1),colspan = 3) 
    plot(times[~flight_mask[:30000]],rwing[~flight_mask[:30000]])
    ax = gca()
    ax.set_ybound((50,90))
    raw_sig = np.mean(dat,axis=1)
    #plot the cycle average of the raw signal
    subplot2grid((n_components+2,5),(1,4))
    cd = cycle_average(raw_sig,cycle_idxs,times)
    plot(cd['time_resamp'],cd['sig_mean'])
    gca().set_yticklabels([])
    #plot an example trace of the raw signal
    subplot2grid((n_components+2,5),(1,1),colspan = 3,sharex = ax) 
    plot(times[~flight_mask[:30000]],raw_sig)
    gca().set_yticklabels([])
    ax = gca()
    #itterate through the components
    for i in range(n_components):
        col = 0 
        row = i
        #first an image of the components
        subplot2grid((n_components+2,5),(row+2,col)) 
        tst_img = np.zeros((256,256))
        tst_img[mask] = A_[:,i] - np.min(A_)
        imshow(tst_img,vmin= 0)#,vmax = 100)
        little_frame = copy.copy(thorax.basis)
        little_frame['p'] = thorax.basis['p']/4
        little_frame['a1'] = thorax.basis['a1']/4
        little_frame['a2'] = thorax.basis['a2']/4
        thorax_view.plot(little_frame,lw = 0.3)
        gca().set_xticks([]); gca().set_yticks([])
        gca().set_xbound((0,256)); gca().set_ybound((0,256))
        #now the cycle average
        subplot2grid((n_components+2,5),(row+2,col+4),colspan = 1)
        cd = cycle_average(sigs[i,:],cycle_idxs,times)
        plot(cd['time_resamp'],cd['sig_mean'])
        gca().set_yticklabels([])
        #now an example trace
        subplot2grid((n_components+2,5),(row+2,col+1),colspan = 3,sharex = ax)
        ax = gca()
        plot(times[:n_samples],sigs[i,:])
        gca().set_yticklabels([])
    ax.set_xbound((40,80))
In [97]:
 
In [97]:
 
In [22]:
 
In [97]:
 
In [82]:
close('all')
In [97]:
figure();plot(ica_sig[1,:]);
In [29]:
 
In [30]:
 
Out[30]:
(3, 4068)
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In [26]:
import numpy as np
from numpy import dot,argsort,diag,where,dstack,zeros, sqrt,float,real,linalg
#from numpy import *
from numpy.linalg import eigh, inv, eig, norm
from numpy.random import randn, rand
import time
try:
    from scipy.stats import skew
    _skew_loaded = True
    from scipy.linalg import orth
    _orth_loaded = True
except:
    _orth_loaded = False
    _skew_loaded = False

## TODOs:
## [ ] show how much variance is accounted by pcs
## [X] function to reshape back to movie
## [X] behave nicely if no scipy available (can't sort by skewness then)
## [ ] simple GUI (traits?)
## [ ] masks from ICs or PCs    

def pca_svd(X):
    ndata, ndim = X.shape
    Y = X - X.mean(axis=-1)[:, np.newaxis] # remove mean
    U,s,Vh = np.linalg.svd(X, full_matrices=False)
    Vh = Vh[:ndata] #only makes sense to return the first num_data
    return U,Vh,s

def pca1 (X, verbose=False):
    """
    Simple principal component analysis
    X as Npix by Nt matrix
    X should be normalized and centered beforehand
    --
    returns:
    - EV (Nt by Nesq:esq>0), matrix of PC 'signals'
     (eigenvalues of temporal covariance matrix). Signals are in columns
    - esq, vector of eigenvalues
    """
    tick = time.clock()
    Nx, Nt = X.shape
    Y = X - X.mean(axis=-1)[:, np.newaxis] # remove mean
    C = dot(Y, Y.T)/Nt # temporal covariance matrix

    es, EV = eigh(C)  # eigenvalues, eigenvectors

    ## take non-negative eigenvalues
    non_neg, = where(es>=0)
    neg = where(es<0)
    if len(neg)>0:
        if verbose:
            print "pca1: Warning, C have %d negative eigenvalues" %len(neg)
        es = es[non_neg]
        EV = EV[:,non_neg]
    #S = diag(np.sqrt(es))
    #whitenmat = dot(EV, inv(S)) # whitenting matrix
    ki = argsort(es)[::-1]      # indices for sorted eigenvalues    
    #U = dot(X, whitenmat)  # spatial filters
    if verbose:
        print "PCA finished in %03.2f sec" %(time.clock() - tick)
    return EV[:,ki], es[ki]

# note: returns EVs stored in columns, e.g.
# [[x1, x2, ... xn],
#  [y1, y2, ... yn],
#  ...]

## Whitening


def whiten_mat1(X):
    EV,es = pca1(X)
    S = diag(np.sqrt(es))
    return dot(EV, inv(S))

def st_ica(X, ncomp = 20,  mu = 0.3):
    """Spatiotemporal ICA for sequences of images
    Input:
    X -- list of 2D arrays or 3D array with first axis = time
    mu -- spatial vs temporal, mu = 0 -> spatial; mu = 1 -> temporal
    """
    data = X
    #data = reshape_from_movie(X) # nframes x npixels
    sh = X[0].shape
    
    pc_filters, pc_signals, ev = whiten_svd(data, ncomp)
    
    mux = sptemp_concat(pc_filters, pc_signals, mu)

    _, W = fastica(mux, whiten=False)
    ica_sig = dot(W, pc_signals)
    ica_filters = dot(dot(diag(1.0/np.sqrt(ev)), W), pc_filters)

    if _skew_loaded:
        skewsorted = argsort(skew(ica_sig, axis=1))[::-1]
    else:
        skewsorted = range(nIC)
    return ica_filters[skewsorted], ica_sig[skewsorted]


def transp(m):
    "conjugate transpose"
    return m.conjugate().transpose()



pow3nonlin = {'g':lambda X: X**3,
              'gprime': lambda X: 3*X**2}

pow3nonlinx = {'g':lambda X,args: X**3,
              'gprime': lambda X,args: 3*X**2}


def _sym_decorrelate(X):
    "W <- W \cdot (W^T \cdot W)^{-1/2}"
    a = dot(X, transp(X))
    ev, EV = linalg.eigh(a)
    
    return dot(dot(dot(EV, np.diag(1.0/np.sqrt(ev))),
                   EV.T),
               X)

def _ica_symm(X, nIC=None, guess=None,
              nonlinfn = pow3nonlin,
              termtol = 5e-7, max_iter = 2e3,
              verbose=False):
    "Simplistic ICA with FastICA algorithm"
    nPC, siglen = map(float, X.shape)
    nIC = nIC or nPC

    guess = guess or np.random.normal(size=(nPC,nPC))
    guess = _sym_decorrelate(guess)

    B, Bprev = zeros(guess.shape, np.float64), guess

    iters, errx = 0,termtol+1
    g, gp = pow3nonlin['g'], pow3nonlin['gprime']

    while (iters < max_iter) and (errx > termtol):
        bdotx = dot(Bprev, X)
        gwtx = g(bdotx)
        gp_wtx = gp(bdotx)/siglen
        B = dot(gwtx, X.T)/siglen - dot(np.diag(gp_wtx.mean(axis=1)), Bprev)
        B = _sym_decorrelate(B)
        errx = max(abs(abs(np.diag(dot(B, Bprev.T)))-1))
        Bprev = np.copy(B)
        iters += 1
    if verbose:
        if iters < max_iter:
            print "Success: ICA Converged in %d tries" %iters
        else:
            print "Fail: reached maximum number of iterations %d reached"%maxiters
    return B.real

def whiten_svd(X, ncomp=None):
    n,p = map(float, X.shape)
    Y = X - X.mean(axis=-1)[:, np.newaxis]
    u, s, vh = linalg.svd(X, full_matrices=False)
    K = (u/s).T[:ncomp]
    return np.dot(K, X), K, s[:ncomp]

def fastica(X, ncomp=None, whiten = True,
            algorithm = 'symmetric',
            tol = 1e-3, max_iter = 1e3, guess = None):
    n,p = map(float, X.shape)
    if whiten:
        XW, K, _ = whiten_svd(X, ncomp) # whitened and projected data
    else:
        XW = X.copy()
    XW *= np.sqrt(p)
    kwargs = {'termtol':tol, 'nonlinfn':pow3nonlin,
              'max_iter':max_iter, 'guess':guess}
    algorithms = {'symmetric':_ica_symm, 'deflation':None}
    fun = algorithms.get(algorithm, 'symmetric')
    W  = fun(XW, **kwargs)
    if whiten:
        S = dot(dot(W,K),X)
    else:
        S = dot(W,X)
    return S, W

def fastica_defl(X, nIC=None, guess=None,
             nonlinfn = pow3nonlin,
             termtol = 5e-7, maxiters = 2e3):
    tick = time.clock()
    nPC, siglen = X.shape
    nIC = nIC or nPC-1
    guess = guess or randn(nPC,nIC)

    if _orth_loaded:
        guess = orth(guess)

    B, Bprev = zeros(guess.shape, np.float64), zeros(guess.shape, np.float64)

    iters, minAbsCos  = 0,0

    errvec = []
    icc = 0
    while icc < nIC:
        w = randn(nPC,1) - 0.5
        w -= dot(dot(B, transp(B)), w)
        w /= norm(w)

        wprev = zeros(w.shape)
        for i in xrange(long(maxiters) +1):
            w -= dot(dot(B, transp(B)), w)
            w /= norm(w)
            #wprev = w.copy()
            if (norm(w-wprev) < termtol) or (norm(w + wprev) < termtol):
                B[:,icc]  = transp(w)
                icc += 1
                break
            wprev = w.copy()
    return B.real, errvec



### Some general utility functions:
### --------------------------------
def gauss_kern(xsize, ysize=None):
    """ Returns a normalized 2D gauss kernel for convolutions """
    xsize = int(xsize)
    ysize = ysize and int(ysize) or xsize
    x, y = mgrid[-xsize:xsize+1, -ysize:ysize+1]
    g = np.exp(-(x**2/float(xsize) + y**2/float(ysize)))
    return g / g.sum()

def gauss_blur(X,size=1.5):
    return signal.convolve2d(X,gauss_kern(size),'same')

def reshape_from_movie(mov):
    l,w,h = mov.shape
    return mov.reshape(l,w*h)

def reshape_to_movie(X,nrows,ncols):
    Np, Nt = X.shape
    return X.T.reshape(Nt,nrows,ncols)

def sptemp_concat(filters, signals, mu):
    if mu == 0:
        out= filters # spatial only
    elif mu == 1:
        out= signals # temporal only
    else:
        out =  np.concatenate(((1-mu)*filters, mu*signals),
                              axis = 1)
    return out / np.sqrt(1-2*mu+2*mu**2)


def DFoF(X):
    """
    Delta F / mean F normalization for each pixel
    """
    xmean = X.mean(axis=1).reshape(-1,1)
    xmean_z = where(xmean==0)
    xmean[xmean_z] = 1.0
    out =  X/xmean - 1
    out[xmean_z,:] = 0
    return out

def DFoSD(X, tslice = None):
    """
    Delta F / S.D.(F) normalization
    """
    if not isinstance(tslice, slice):
        tslice = tslice and slice(*tslice) or slice(None)
    xsd = X[:,tslice].std(axis=1).reshape(-1,1)
    xmean = X[:,tslice].mean(axis=1).reshape(-1,1)
    return (X-xmean)/xsd

def demean(X):
    """
    Remove mean over time from each pixel
    Frames are flattened and are in columns
    """
    xmean = X.mean(axis=0).reshape(1,-1)
    return X - xmean


def winvhalf(X):
    "raise matrix to power -1/2"
    e, V = eigh(X)
    return dot(V, dot(diag((e+0j)**-0.5),transp(V)))

def msqrt(X):
    e, V = eigh(X)
    return dot(V, dot(diag(((e+0j)**0.5)), transp(V)))

def mpower(M, p):
    """
    Matrix exponentiation, works for Hermitian matrices
    """
    e,EV = linalg.eigh(M)
    return dot(transp(EV),
               dot(diag((e+0j)**p), EV))


def mask4overlay(mask,colorind=0):
    """
    Put a binary mask in some color channel
    and make regions where the mask is False transparent
    """
    sh = mask.shape
    z = np.zeros(sh)
    stack = dstack((z,z,z,np.ones(sh)*mask))
    stack[:,:,colorind] = mask
    return stack


def flcompose2(f1,f2):
    "Compose two functions from left to right"
    def _(*args,**kwargs):
        return f2(f1(*args,**kwargs))
    return _
                  
def flcompose(*funcs):
    "Compose a list of functions from left to right"
    return reduce(flcompose2, funcs)
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In []:
 
In [184]:
nc = 3

#ica_data = warped_movie[~flight_mask[:10000],:,:]

S_,A_,C_ = runICA(ica_data,b3_mask,n_components = nc)

plot_ICA(warped_movie,
        b3_mask,
        flight_mask,
        times,
        rwing,
        cycle_idxs,
        S_,A_,C_)

plot_comp_img(A_,b3_mask)
In [190]:
print mean(A_,axis = 0)
[ 10.70822287 -87.15541347  79.43696963]

In [165]:
 
In [165]:
 
In [165]:
 
In [165]:
 
In []:
 
In [19]:
nc = 12
all_muscle_mask = b1_mask | b2_mask \
              | b3_mask | i1_mask \
              | i2_mask | iii1_mask \
              | iii3_mask | iii24_mask \
              | hg1_mask | hg2_mask \
              | hg3_mask | hg4_mask

ica_data = warped_movie[~flight_mask[:5000],:,:]

S_,A_,C_ = runICA(ica_data,all_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
        all_muscle_mask,
        flight_mask,
        times,
        rwing,
        cycle_idxs,
        S_,A_,C_)
plot_comp_img(A_,all_muscle_mask)
/home/flyranch/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/sklearn/decomposition/fastica_.py:110: UserWarning: FastICA did not converge. You might want to increase the number of iterations.
  ' to increase the number of iterations.')

In [31]:
nc = 4
anterior_muscle_mask = b1_mask | b2_mask \
              | b3_mask | i1_mask \

ica_data = warped_movie[~flight_mask[:5000],:,:]

S_,A_,C_ = runICA(ica_data,anterior_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
        anterior_muscle_mask,
        flight_mask,
        times,
        rwing,
        cycle_idxs,
        S_,A_,C_)

figure()
plot_comp_img(A_,anterior_muscle_mask)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-31-0fb8a4dcebf1> in <module>()
     13 
     14 figure()
---> 15 plot_comp_img(A_,anterior_muscle_mask)

<ipython-input-29-ab0b08239d63> in plot_comp_img(A_, mask)
     17         col = mod(i,3)
     18         row = i/3
---> 19         subplot2grid((n_components/3,3),(row,col))
     20         tst_img = np.zeros((256,256))
     21         tst_img[mask] = A_[:,i] - np.min(A_)

/home/flyranch/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/pyplot.pyc in subplot2grid(shape, loc, rowspan, colspan, **kwargs)
   1171                                                    rowspan=rowspan,
   1172                                                    colspan=colspan)
-> 1173     a = fig.add_subplot(subplotspec, **kwargs)
   1174     bbox = a.bbox
   1175     byebye = []

/home/flyranch/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/figure.pyc in add_subplot(self, *args, **kwargs)
    956                     self._axstack.remove(ax)
    957 
--> 958             a = subplot_class_factory(projection_class)(self, *args, **kwargs)
    959 
    960         self._axstack.add(key, a)

/home/flyranch/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/axes/_subplots.pyc in __init__(self, fig, *args, **kwargs)
     73             raise ValueError('Illegal argument(s) to subplot: %s' % (args,))
     74 
---> 75         self.update_params()
     76 
     77         # _axes_class is set in the subplot_class_factory

/home/flyranch/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/axes/_subplots.pyc in update_params(self)
    113         self.figbox, self.rowNum, self.colNum, self.numRows, self.numCols = \
    114             self.get_subplotspec().get_position(self.figure,
--> 115                                                 return_all=True)
    116 
    117     def is_first_col(self):

/home/flyranch/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/gridspec.pyc in get_position(self, fig, return_all)
    427 
    428         rowNum, colNum =  divmod(self.num1, ncols)
--> 429         figBottom = figBottoms[rowNum]
    430         figTop = figTops[rowNum]
    431         figLeft = figLefts[colNum]

IndexError: list index out of range
In [32]:
nc = 2

ica_data = warped_movie[~flight_mask[:5000],:,:]
b3_bool_mask = b3_mask & ~(i1_mask | b1_mask | b2_mask)

S_,A_,C_ = runICA(warped_movie,b3_bool_mask,n_components = 2)
plot_ICA(warped_movie,
        b3_bool_mask,
        flight_mask,
        times,
        rwing,
        cycle_idxs,
        S_,A_,C_)
#display(gcf());close()
In [44]:
fly = U_X_J.flies[0]
#load the warped movie
warped_movie = tifffile.TiffFile(fly.fly_path + 'warped_imgs.tiff').asarray()
#only did one experiment on these flies
expmnt = fly.experiments.values()[0]
#Map some variables
stim_cond = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['DAC2'])
lwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph0'])
rwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])
lmr = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph2'])
aux = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph3'])
xpos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Xpos'])
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])

flight_mask = aux<0.0
trial_selection_mask = stim_cond>0
cycle_idxs = collect_cycles(ypos,stim_cond,trial_selection_mask)
In [45]:
nc = 4
anterior_muscle_mask = b1_mask | b2_mask \
              | b3_mask | i1_mask \

ica_data = warped_movie[~flight_mask[:5000],:,:]

S_,A_,C_ = runICA(ica_data,anterior_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
        anterior_muscle_mask,
        flight_mask,
        times,
        rwing,
        cycle_idxs,
        S_,A_,C_)

figure()
#plot_comp_img(A_,anterior_muscle_mask)
Out[45]:
<matplotlib.figure.Figure at 0xf41b290>
In [46]:
nc = 3
#i1_mask = b1_mask | b2_mask \
#              | b3_mask | i1_mask \

ica_data = warped_movie[~flight_mask[:5000],:,:]

S_,A_,C_ = runICA(ica_data,i1_mask,n_components = nc)
plot_ICA(warped_movie,
        i1_mask,
        flight_mask,
        times,
        rwing,
        cycle_idxs,
        S_,A_,C_)

figure()
#plot_comp_img(A_,anterior_muscle_mask)
Out[46]:
<matplotlib.figure.Figure at 0x544ef290>
In [48]:
nc = 3
#i1_mask = b1_mask | b2_mask \
#              | b3_mask | i1_mask \

ica_data = warped_movie[~flight_mask[:5000],:,:]

S_,A_,C_ = runICA(ica_data,b2_mask,n_components = nc)
plot_ICA(warped_movie,
        b2_mask,
        flight_mask,
        times,
        rwing,
        cycle_idxs,
        S_,A_,C_)

figure()
#plot_comp_img(A_,anterior_muscle_mask)
Out[48]:
<matplotlib.figure.Figure at 0x881ffd0>
In [47]:
draw()
In [34]:
plot(proj[1,:])
Out[34]:
[<matplotlib.lines.Line2D at 0x2f7297d0>]
In [57]:
plot(cycle_average(S_[:,0],cycle_idxs,times)['sig_median'])
Out[57]:
[<matplotlib.lines.Line2D at 0x150e36350>]
In [126]:
s1 = sum(warped_movie[:,mask],axis = 1)
In [122]:
i1_bool_mask = i1_mask & ~(b2_mask)
mask = cv2.pyrDown(cv2.pyrDown(np.float64(i1_bool_mask))) >0
s2 = sum(warped_movie[:,mask],axis = 1)
In [123]:
from sklearn.decomposition import FastICA
ica = FastICA(n_components=1)
S  = c_[s1,s2]/c_[s1,s2].std(axis = 0)
S_ = ica.fit_transform(c_[s2,s1])  # Reconstruct signals
A_ = ica.mixing_  # Get estimated mixing matrix
C_ = ica.components_
In [127]:
imshow(b3_bool_mask)
Out[127]:
<matplotlib.image.AxesImage at 0x173e5fb50>
In [119]:
plot(flight_mask)
Out[119]:
[<matplotlib.lines.Line2D at 0x150b5ea10>]

This approach seems like it is useful in the sense that it identifies neworks of activity, but it is difficult to interpret in terms of isolating the signal from individual muscles. This is likely due to the fact that optimization of the mixing matrix is not constrained by any spatial information. I have found several papers that describe ICA with spatial constraints, but implementing these require that I roll my own ICA optimization algorithm - and this would take some time. As a first pass I took an ad-hoc approach to applying spatial constraints: first by breaking the analysis up into the anterior and posterior cluster of muscles, then looking at the components of each muscle one-by-one.

In [15]:
#the anterior cluster
nc = 4
anterior_muscle_mask = b1_mask | b2_mask \
                       |b3_mask | i1_mask
S_,A_ = runICA(anterior_muscle_mask,n_components = nc)
plot_ICA(anterior_muscle_mask,S_,A_,n_components = nc);display(gcf());close()
In [17]:
nc = 3
b2_bool_mask =  b2_mask 
S_,A_ = runICA(b2_bool_mask,n_components = nc)
plot_ICA(b2_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [18]:
nc = 3
i1_bool_mask =  i1_mask
S_,A_ = runICA(i1_bool_mask,n_components = nc)
plot_ICA(i1_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [19]:
nc = 3
b3_bool_mask =  b3_mask
S_,A_ = runICA(b3_bool_mask,n_components = nc)
plot_ICA(b3_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [20]:
nc = 3
b1_bool_mask =  b1_mask
S_,A_ = runICA(b1_bool_mask,n_components = nc)
plot_ICA(b1_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [21]:
nc = 3
iii3_bool_mask =  iii3_mask 
S_,A_ = runICA(iii3_bool_mask,n_components = nc)
plot_ICA(iii3_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [22]:
nc = 3
iii24_bool_mask =  iii24_mask 
S_,A_ = runICA(iii24_bool_mask,n_components = nc)
plot_ICA(iii24_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [23]:
nc = 3
hg1_bool_mask =  hg1_mask 
S_,A_ = runICA(hg1_bool_mask,n_components = nc)
plot_ICA(hg1_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [24]:
nc = 3
hg2_bool_mask =  hg2_mask 
S_,A_ = runICA(hg2_bool_mask,n_components = nc)
plot_ICA(hg2_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [25]:
nc = 3
hg3_bool_mask =  hg3_mask 
S_,A_ = runICA(hg3_bool_mask,n_components = nc)
plot_ICA(hg3_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [26]:
nc = 3
hg4_bool_mask =  hg4_mask 
S_,A_ = runICA(hg4_bool_mask,n_components = nc)
plot_ICA(hg4_bool_mask,S_,A_,n_components = nc);display(gcf());close()

it is also possible to define masks that atempt to excude pixels from overlaping regions

In [27]:
nc = 3
b2_only_boolmask =  i1_mask & ~(b2_mask | b1_mask |b3_mask)
S_,A_ = runICA(b2_only_boolmask,n_components = nc)
plot_ICA(b2_only_boolmask,S_,A_,n_components = nc);display(gcf());close()
In [0]:
 
In [0]:
 
In [64]:
 
In [45]:
### function to itterate through a large list of flies
### and make figures of the fits - probably should migrate this to
### an independent script
import db_access as dba
import flylib as flb
import cv2
import cPickle
fly_db = dba.get_db()

def plot_maxfits(flylist):
    numfigs = np.ceil(len(flylist)/6.)
    chunks = np.array_split(flylist,numfigs)
    for chunk in chunks:
        fit_swarm = flb.Squadron(fly_db,chunk)
        figure(figsize = (8.5,11))
        for i,fly in enumerate(fit_swarm.flies):
            output_shape = (1024,1024)
            #fly = swarm.flies[0]
            tiffdata = np.array(fly.fly_record['experiments'].values()[0]['tiff_data']['images'])
            pkname = fly.fly_path + '/basis_fits.cpkl'
            f = open(pkname)
            img_data = cPickle.load(f)
            f.close()
            test_basis = Basis()
            for key in ['A','p','a1','a2']:
                test_basis[key] = img_data[key]

            src_p0 = test_basis['a1'] + test_basis['p']
            src_p1 = test_basis['p']
            src_p2 = test_basis['a2'] + test_basis['p']
            dst_p0 = confocal_basis['a1'] + confocal_basis['p']
            dst_p1 = confocal_basis['p']
            dst_p2 = confocal_basis['a2'] + confocal_basis['p']
            maximg = np.max(tiffdata,axis=0)
            col = mod(i,2)
            row = i/2
            subplot2grid((3,2),(row,col))
            A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
            imshow(cv2.warpAffine(maximg.T,A,output_shape),cmap = cm.gray)
            thorax_view.plot(thorax.basis,alpha = 0.5)
            gca().set_ybound(0,1024)
            gca().set_xbound(0,1024)
            gca().set_xticks([])
            gca().set_yticks([])
            gca().set_title(('Fly%04d')%(fly.fly_num))